package twoD.IsoPfem;

import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.ELAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;
import inf.jlinalg.Vector2D;
import twoD.IsoPfem.Structure.CType;
import twoD.IsoPfem.math.FMath;
import twoD.IsoPfem.math.FunctionRToR;
import twoD.IsoPfem.math.Integrator;
import twoD.IsoPfem.math.Interval;

/**
 * implementation class for linear triangle
 * 
 * @author hbui
 * 
 */
public class CST extends Element2D {

	public CST(IGeometry T) {
		super(T);
		setCType(CType.PLANE_STRESS);
	}

	public CST(RealTable realTable, IGeometry T) {
		this(T);
		this.setRealTable(realTable);
	}

	@Override
	public double getVolume() {
		Vector2D p1 = getNodes(0).getPosition();
		Vector2D p2 = getNodes(1).getPosition();
		Vector2D p3 = getNodes(2).getPosition();
		Vector2D d13 = p1.subtract(p3);
		Vector2D d23 = p2.subtract(p3);
		double dx1n13 = d13.getX1();
		double dx1n23 = d23.getX1();
		double dx2n13 = d13.getX2();
		double dx2n23 = d23.getX2();
		double detJ = dx1n13 * dx2n23 - dx1n23 * dx2n13;

		return getThickness() * Math.abs(detJ) / 2;
	}

	public IMatrix getBMatrix() {
		Array2DMatrix b = new Array2DMatrix(3, 6);
		Vector2D p1 = getNodes(0).getPosition();
		Vector2D p2 = getNodes(1).getPosition();
		Vector2D p3 = getNodes(2).getPosition();
		Vector2D d12 = p1.subtract(p2);
		Vector2D d13 = p1.subtract(p3);
		Vector2D d23 = p2.subtract(p3);
		double dx1n12 = d12.getX1();
		double dx1n13 = d13.getX1();
		double dx1n23 = d23.getX1();
		double dx2n12 = d12.getX2();
		double dx2n13 = d13.getX2();
		double dx2n23 = d23.getX2();
		double detJ = dx1n13 * dx2n23 - dx1n23 * dx2n13;

		// b-matrix
		b.set(0, 0, dx2n23 / detJ);
		b.set(0, 2, -dx2n13 / detJ);
		b.set(0, 4, dx2n12 / detJ);
		b.set(1, 1, -dx1n23 / detJ);
		b.set(1, 3, dx1n13 / detJ);
		b.set(1, 5, -dx1n12 / detJ);
		b.set(2, 0, -dx1n23 / detJ);
		b.set(2, 1, dx2n23 / detJ);
		b.set(2, 2, dx1n13 / detJ);
		b.set(2, 3, -dx2n13 / detJ);
		b.set(2, 4, -dx1n12 / detJ);
		b.set(2, 5, dx2n12 / detJ);
		return b;
	}

	@Override
	public IVector computeInternalForceVector() {
		// IMatrix ke = computeTangentStiffnessMatrix();
		// IVector ue = getUe();
		// IVector rie = new ArrayVector(ue.getSize());
		// BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, ke, ue, 0.0, rie);

		// using the method below will ensure derived type of element (i.e
		// crack) can performed correctly since this is the general equation for
		// sigma calculation for constant strain triangle
		IVector rie = new ArrayVector(6);
		IMatrix B = getBMatrix();
		IVector sigma = computeSigma();
		BLAM.multiply(getVolume(), BLAM.TRANSPOSE, B, sigma, 0.0, rie);
		return rie;
	}

	@Override
	public IMatrix computeTangentStiffnessMatrix() {
		IMatrix ke = new Array2DMatrix(6, 6);
		double A = getVolume();
		ELAM.btCB(getBMatrix(), getCMatrix(), A, ke);
		// System.out.println(this.getId() + ":"
		// + MatrixFormat.format(getBMatrix()));
		// System.out
		// .println("ke " + this.getId() + "=" + MatrixFormat.format(ke));
		return ke;
	}

	@Override
	public void computeExternalForceVector() {
		Vector2D p1 = getNodes(0).getPosition();
		Vector2D p2 = getNodes(1).getPosition();
		Vector2D p3 = getNodes(2).getPosition();
		Vector2D d13 = p1.subtract(p3);
		Vector2D d23 = p2.subtract(p3);
		double dx1n13 = d13.getX1();
		double dx1n23 = d23.getX1();
		double dx2n13 = d13.getX2();
		double dx2n23 = d23.getX2();
		double absDetJt = getThickness()
				* Math.abs(dx1n13 * dx2n23 - dx1n23 * dx2n13);

		FunctionRToR N1 = new FunctionRToR() {

			@Override
			public double valueAt(double x) {
				return 1 - x;
			}

			@Override
			public int getP() {
				return 1;
			}

			@Override
			public Interval getDomain() {
				return new Interval(0, 1);
			}

			@Override
			public double[] getBreakpoints() {
				return null;
			}

			@Override
			public double derivativeAt(double x) {
				return -1;
			}
		};
		FunctionRToR N2 = new FunctionRToR() {

			@Override
			public double valueAt(double x) {
				return x;
			}

			@Override
			public int getP() {
				return 1;
			}

			@Override
			public Interval getDomain() {
				return new Interval(0, 1);
			}

			@Override
			public double[] getBreakpoints() {
				return null;
			}

			@Override
			public double derivativeAt(double x) {
				return 1;
			}
		};

		for (EdgeLoad load : edgeLoad_) {
			Node n1 = load.getN1();
			Node n2 = load.getN2();
			double f1 = Integrator
					.integrate(FMath.multiply(N1, load.getLoad()));
			double f2 = Integrator
					.integrate(FMath.multiply(N2, load.getLoad()));
			Vector2D dn = n2.getPosition().subtract(n1.getPosition());
			double l = dn.norm2();
			double c = dn.getX1() / l;
			double s = dn.getX2() / l;
			/*
			 * update nodal force in n1 and n2
			 */
			// n1.setNodalForce(new Force(-f1 * s * absDetJt, f1 * c *
			// absDetJt));
			// n2.setNodalForce(new Force(-f2 * s * absDetJt, f2 * c *
			// absDetJt));
			n1.addNodalForce(-f1 * s * absDetJt, f1 * c * absDetJt);
			n2.addNodalForce(-f2 * s * absDetJt, f2 * c * absDetJt);
		}
	}

	@Override
	public IVector computeEpsilon(double... xi) {
		IMatrix B = getBMatrix();
		IVector Ue = this.getUe();
		// System.out
		// .println("Ue" + this.getId() + " =" + MatrixFormat.format(Ue));
		IVector e = new ArrayVector(3);
		BLAM.multiply(1.0, BLAM.NO_TRANSPOSE, B, Ue, 0.0, e);
		return e;
	}

}
